clear
current_directory = pwd;

num_comps = 16
parpool('local',num_comps)

rng(17)

% Load productivity shocks

avg_literature_ag_shocks = importdata('avg_literature_ag_shocks.mat');
Cline_shocks = avg_literature_ag_shocks;
manu_shocks = importdata('manu_productivity_shocks_2080_2099_CSIRO_climateadaptation.mat')';
manu_shocks = ones(1,length(manu_shocks))+manu_shocks;
serv_shocks = importdata('serv_productivity_shocks_2080_2099_CSIRO_climateadaptation.mat')';
serv_shocks = ones(1,length(serv_shocks))+serv_shocks;

% First, run the baseline simulation
baselinesimulation

% Save all the moments from the baseline simulation
% Import shares and domestic expenditure shares for all countries and bilateral pairs
ag_import_shares_mat_baseline = ag_import_share_mat;
ag_domestic_expenditures_vec_baseline = ag_import_share_mat(logical(eye(n_countries)));
manu_import_shares_mat_baseline = manu_import_share_mat;
manu_domestic_expenditures_vec_baseline = manu_import_share_mat(logical(eye(n_countries)));
ag_trade_mat_baseline = ag_trade_mat;
manu_trade_mat_baseline = manu_trade_mat;

% GDP shares by sector 
ag_GDP_share_baseline = sim_ag_shares_GDP;
manu_GDP_share_baseline = sim_manu_shares_GDP;
serv_GDP_share_baseline = sim_serv_shares_GDP;
ag_exp_share_baseline = exp_share_ag_vec;
manu_exp_share_baseline = exp_share_manu_vec;
serv_exp_share_baseline = exp_share_serv_vec;
ag_net_exports_baseline = ag_GDP_share_baseline-ag_exp_share_baseline;

% Aggregate sectoral consumption and price indices by country
C_vec_baseline = C_vec;
P_vec_baseline = P_vec;
P_a_vec_baseline = P_a_vec;
P_m_vec_baseline = P_m_vec;
P_s_vec_baseline = P_s_vec;
C_a_vec_baseline = C_a_vec;
C_m_vec_baseline = C_m_vec;
C_s_vec_baseline = C_s_vec;
income_per_capita_baseline = wage_guesses;
total_income_vec_baseline = total_income_vec;
GDP_baseline = P_a_vec_baseline.*C_a_vec_baseline+P_m_vec_baseline.*C_m_vec_baseline+P_s_vec_baseline.*C_s_vec_baseline;
P_Tornqvist_vec_baseline = P_a_vec_baseline.^(ag_exp_share_baseline).*P_m_vec_baseline.^(manu_exp_share_baseline)...
    .*P_s_vec_baseline.^(serv_exp_share_baseline);
GDP_Tornqvist_baseline = total_income_vec_baseline./P_Tornqvist_vec_baseline;

% Choose trade policy counterfactual permutation - no tariffs, free trade
% agreement, no red tape barriers, and removing all
for trade_counterfactual = 1:1:3
    
if trade_counterfactual == 1
ag_trade_cost_mat = csvread('ag_trade_costs_no_tariffs_FTA.csv', 1, 0);
manu_trade_cost_mat = csvread('manu_trade_costs_no_tariffs_FTA.csv', 1, 0);
T_multiplier = importdata('temp_T_multiplier_lowtradecost_no_tariffs_FTA.mat');
elseif trade_counterfactual == 2
ag_trade_cost_mat = csvread('ag_trade_costs_no_red_tape.csv', 1, 0);
manu_trade_cost_mat = csvread('manu_trade_costs_no_red_tape.csv', 1, 0);
T_multiplier = importdata('temp_T_multiplier_lowtradecost_no_red_tape.mat');
elseif trade_counterfactual == 3
ag_trade_cost_mat = csvread('ag_trade_costs_no_policy_barriers.csv', 1, 0);
manu_trade_cost_mat = csvread('manu_trade_costs_no_policy_barriers.csv', 1, 0);
T_multiplier = importdata('temp_T_multiplier_lowtradecost_no_policy_barriers.mat');
end

baselinesimulation_lowtradecost_policy_counterfactuals

% Save all the same moments as above, but label them as the low trade cost scenario, without climate change

ag_import_shares_mat_baseline_lowtradecost = ag_import_share_mat;
ag_domestic_expenditures_vec_baseline_lowtradecost = ag_import_share_mat(logical(eye(n_countries)));
manu_import_shares_mat_baseline_lowtradecost = manu_import_share_mat;
manu_domestic_expenditures_vec_baseline_lowtradecost = manu_import_share_mat(logical(eye(n_countries)));
ag_trade_mat_baseline_lowtradecost = ag_trade_mat;
manu_trade_mat_baseline_lowtradecost = manu_trade_mat;
ag_GDP_share_baseline_lowtradecost = sim_ag_shares_GDP;
manu_GDP_share_baseline_lowtradecost = sim_manu_shares_GDP;
serv_GDP_share_baseline_lowtradecost = sim_serv_shares_GDP;
ag_exp_share_baseline_lowtradecost = exp_share_ag_vec;
manu_exp_share_baseline_lowtradecost = exp_share_manu_vec;
serv_exp_share_baseline_lowtradecost = exp_share_serv_vec;
ag_net_exports_baseline_lowtradecost = ag_GDP_share_baseline_lowtradecost-ag_exp_share_baseline_lowtradecost;
C_vec_baseline_lowtradecost = C_vec;
P_vec_baseline_lowtradecost = P_vec;
P_a_vec_baseline_lowtradecost = P_a_vec;
P_m_vec_baseline_lowtradecost = P_m_vec;
P_s_vec_baseline_lowtradecost = P_s_vec;
C_a_vec_baseline_lowtradecost = C_a_vec;
C_m_vec_baseline_lowtradecost = C_m_vec;
C_s_vec_baseline_lowtradecost = C_s_vec;
income_per_capita_baseline_lowtradecost = wage_guesses;
total_income_vec_baseline_lowtradecost = total_income_vec;
GDP_baseline_lowtradecost = P_a_vec_baseline_lowtradecost.*C_a_vec_baseline_lowtradecost+P_m_vec_baseline_lowtradecost.*C_m_vec_baseline_lowtradecost+P_s_vec_baseline_lowtradecost.*C_s_vec_baseline_lowtradecost;
P_Tornqvist_vec_baseline_lowtradecost = P_a_vec_baseline_lowtradecost.^(ag_exp_share_baseline_lowtradecost).*P_m_vec_baseline_lowtradecost.^(manu_exp_share_baseline_lowtradecost)...
    .*P_s_vec_baseline_lowtradecost.^(serv_exp_share_baseline_lowtradecost);
GDP_Tornqvist_baseline_lowtradecost = total_income_vec_baseline_lowtradecost./P_Tornqvist_vec_baseline_lowtradecost;


% Now simulate the low trade cost scenario *with* climate change

ag_T_vec = ag_T_vec.*Cline_shocks;
manu_T_vec = manu_T_vec.*manu_shocks;
serv_T_vec = serv_T_vec.*serv_shocks;

sim_Frechets
P_a_vec = geomean(ag_consumed_prices_mat);
P_m_vec = geomean(manu_consumed_prices_mat);
P_s_vec = geomean(serv_consumed_prices_mat);

find_equilibrium_parallel_lowtradecost

sim_ag_shares_GDP = ag_income_vec./total_income_vec;
sim_manu_shares_GDP = manu_income_vec./total_income_vec;
sim_serv_shares_GDP = serv_income_vec./total_income_vec;

% Again save all the moments, but this time for the low trade cost scenario with climate change

ag_import_shares_mat_counterfactual_lowtradecost = ag_import_share_mat;
ag_domestic_expenditures_vec_counterfactual_lowtradecost = ag_import_share_mat(logical(eye(n_countries)));
manu_import_shares_mat_counterfactual_lowtradecost = manu_import_share_mat;
manu_domestic_expenditures_vec_counterfactual_lowtradecost = manu_import_share_mat(logical(eye(n_countries)));
ag_trade_mat_counterfactual_lowtradecost = ag_trade_mat;
manu_trade_mat_counterfactual_lowtradecost = manu_trade_mat;
ag_GDP_share_counterfactual_lowtradecost = sim_ag_shares_GDP;
manu_GDP_share_counterfactual_lowtradecost = sim_manu_shares_GDP;
serv_GDP_share_counterfactual_lowtradecost = sim_serv_shares_GDP;
ag_exp_share_counterfactual_lowtradecost = exp_share_ag_vec;
manu_exp_share_counterfactual_lowtradecost = exp_share_manu_vec;
serv_exp_share_counterfactual_lowtradecost = exp_share_serv_vec;
ag_net_exports_counterfactual_lowtradecost = ag_GDP_share_counterfactual_lowtradecost-ag_exp_share_counterfactual_lowtradecost;
C_vec_counterfactual_lowtradecost = C_vec;
P_vec_counterfactual_lowtradecost = P_vec;
P_a_vec_counterfactual_lowtradecost = P_a_vec;
P_m_vec_counterfactual_lowtradecost = P_m_vec;
P_s_vec_counterfactual_lowtradecost = P_s_vec;
C_a_vec_counterfactual_lowtradecost = C_a_vec;
C_m_vec_counterfactual_lowtradecost = C_m_vec;
C_s_vec_counterfactual_lowtradecost = C_s_vec;
income_per_capita_counterfactual_lowtradecost = wage_guesses;
total_income_vec_counterfactual_lowtradecost = total_income_vec;
GDP_counterfactual_lowtradecost = P_a_vec_counterfactual_lowtradecost.*C_a_vec_counterfactual_lowtradecost+P_m_vec_counterfactual_lowtradecost.*C_m_vec_counterfactual_lowtradecost+P_s_vec_counterfactual_lowtradecost.*C_s_vec_counterfactual_lowtradecost;
% P_Tornqvist_vec_counterfactual_lowtradecost = P_a_vec_counterfactual_lowtradecost.^(ag_exp_share_counterfactual_lowtradecost).*P_m_vec_counterfactual_lowtradecost.^(manu_exp_share_counterfactual_lowtradecost)...
%     .*P_s_vec_counterfactual_lowtradecost.^(serv_exp_share_counterfactual_lowtradecost);
P_Tornqvist_vec_counterfactual_lowtradecost = P_a_vec_counterfactual_lowtradecost.^((ag_exp_share_baseline_lowtradecost+ag_exp_share_counterfactual_lowtradecost)/2)...
    .*P_m_vec_counterfactual_lowtradecost.^((manu_exp_share_baseline_lowtradecost+manu_exp_share_counterfactual_lowtradecost)/2)...
    .*P_s_vec_counterfactual_lowtradecost.^((serv_exp_share_baseline_lowtradecost+serv_exp_share_counterfactual_lowtradecost)/2);
GDP_Tornqvist_counterfactual_lowtradecost = total_income_vec_counterfactual_lowtradecost./P_Tornqvist_vec_counterfactual_lowtradecost;


% Run a function to calculate equivalent variation in the low trade cost scenario 

EV_lowtradecost = calc_EV(C_vec_counterfactual_lowtradecost, income_per_capita_baseline_lowtradecost, P_a_vec_baseline_lowtradecost, P_m_vec_baseline_lowtradecost, P_s_vec_baseline_lowtradecost, ...
    epsilon_a,epsilon_m,epsilon_s,sigma,omega_a,omega_m,omega_s,num_comps,n_countries);

% Save all the matrices and export data as results

GDP_Tornqvist_mat = [GDP_Tornqvist_baseline' GDP_Tornqvist_baseline_lowtradecost' GDP_Tornqvist_counterfactual_lowtradecost'];
GDP_Tornqvist_changes_mat = [GDP_Tornqvist_baseline_lowtradecost'./GDP_Tornqvist_baseline' ...
    GDP_Tornqvist_counterfactual_lowtradecost'./GDP_Tornqvist_baseline_lowtradecost'];

ag_GDP_share_mat = [ag_GDP_share_baseline' ag_GDP_share_baseline_lowtradecost' ag_GDP_share_counterfactual_lowtradecost'];

ag_domestic_expenditures_mat = [ag_domestic_expenditures_vec_baseline ag_domestic_expenditures_vec_baseline_lowtradecost ag_domestic_expenditures_vec_counterfactual_lowtradecost];

food_prices_mat = [P_a_vec_baseline' P_a_vec_baseline_lowtradecost' P_a_vec_counterfactual_lowtradecost'];

outputs = [GDP_Tornqvist_changes_mat EV_lowtradecost' ag_GDP_share_mat ag_domestic_expenditures_mat  ...
    ag_net_exports_baseline' ag_net_exports_baseline_lowtradecost' ag_net_exports_counterfactual_lowtradecost' food_prices_mat];


if trade_counterfactual == 1
dlmwrite('results_lowtradecost_no_tariffs_FTA.csv',outputs)
elseif trade_counterfactual == 2
dlmwrite('results_lowtradecost_no_red_tape.csv',outputs)
elseif trade_counterfactual == 3
dlmwrite('results_lowtradecost_no_policy_barriers.csv',outputs)
end

end
